The present invention relates to a method for automatically generating a numerical calculation program by a finite difference method for a physical model which is expressed by an arbitrary differential equation, and more particularly, to a program generating method suitable for fluid analysis and a region dividing method therefor.
With development of digital computers, there have been proposed a variety of methods for approximately solving numerical calculations in a practically acceptable accuracy. All of these methods are based on solving a finite difference equation which is a substitution of a partial differential equation with a finite difference, and referred to as a finite difference method. According to these methods, a spatial region is fractionated into a mesh form (generally a square mesh form) and a point in the spatial region is represented by a mesh point. In numerical calculation methods for partial differential equations, when a differential equation is established in an object region for the simulation, the object region is divided into small parts in a mesh form, wherein the partial differential equation is discretized at the mesh point, and a numerical solution is obtained by solving this discrete equation. The discrete equation, mentioned above, refers to a relational equation of a variable at a discrete point at which a differential operator can be expanded, and the expansion of the differential operator is called discretization.
On the other hand, a method for generating a program by encoding a discrete equation for executing the finite difference method is described, for example, in "Numerical Simulation Terms DEQSOL" by Sagawa et al in Transactions of National Conference of Information Processing Society of Japan, "Computer System" Symposium, Vol. 30, No. 1, January 1989. This method generates a program by introducing discrete equations in accordance with a discretization formula derived by Taylor expansion from partial differential equations at inner meshes as well as from boundary conditions on boundary meshes, and then encoding those discrete equations.
Conventionally, as methods for automatically generating programs by the finite difference method, there have been proposed, for example, an automatic numerical calculation method described in U.S. Pat. No. 4,841,479 assigned to the present assignee and issued Jun. 20, 1989 (corresponding to JA-A-62-212771), and "Program Generating Method" described in U.S. Ser. No. 361,283 (corresponding to JP-A-1-307826), assigned to the present assignee and filed Jun. 5, 1989. There have also been published papers such as "The Computation of Compressible and Incompressible Recirculating Flows by a Non-iterative Implicit Scheme" by R. I. Issa in Journal of Computer Physics, 62, 66-82 (1986) and "Comparison of Third Order Upwind-Difference Methods in High Reinolds Number" by Kawai et al in Transactions of the Japan Society of Mechanical Engineers, B-52-477 (May 1986).
According to the specification and drawings of the above referred U.S. Pat. No. 4,841,479 (U.S. Ser. No. 16406), the method for automatically generating a program by the finite difference method automatically generates a calculation program such that, assuming that a partial differential equation, an object region or domain for the simulation, mesh information, a variable and a constant are given as inputted information, evaluation points (discrete points) for the variable and the constant are located at mesh points, and a discretization satisfying the conservation law of a physical quantity is performed by control volumes surrounding the discrete points. The control volume is provided for dividing the object region for the simulation into minute regions, the centers of which are located at the respective coordinate points of the discrete points, and performing processing by using these minute regions. Stated other way, the thinking of the control volume is based on that a quantity inputted to a mesh must be equal to a quantity outputted from the mesh, and the value at each discrete point is considered to be a value derived by an integration in the minute region and also a value that is finally divided by the dimension of the minute region. According to this thinking, the total sum of the values obtained in these minute regions becomes equal to the value on the whole object region. It has been known, as shown in the foregoing document "Journal of Computer PHYSICS" and U.S. Ser. No. 361,283, that the discretization by using the control volume exhibits a higher accuracy than a discretization by using Taylor expansion. However, the above method cannot generate simulation programs relative to fluid problems due to mathematical characteristics of a field equation. This is because the above method does not include respective numerical calculation methods for a staggered mesh and an upwind difference necessary for computational fluid dynamics.
The staggered mesh in the above context is a method for deriving a discrete equation by locating a discrete point of a variable at a position different from a mesh point, and this is regarded as the standard since it was published in 1964. It can be also said that the staggered mesh is a method for positioning mesh points which is effective in solving fluid problems dominated by the Navier-Stokes equation in a numerical analytical manner, as shown in the above-mentioned document "Journal of Computer PHYSICS".
The upwind difference refers to a method relative to discretization of a differential operator, wherein a proportion of discrete values on the upstream side of a fluid is chosen to be a larger value while a proportion of discrete values on the downstream side of the fluid is chosen to be a smaller value, and also this is a technique for more stably obtaining a solution in a numerical manner for analytically solving a fluid problem.
As for the upwind difference, as described in the foregoing document "Comparison of Third Order Upwind-Difference Methods in High Reinolds Number" by Kawai et al in Transactions of the Japan Society of Mechanical Engineers, B-52-477 (May 1986), there are a variety of upwind methods. The upwind difference is also regarded as the standard in the computational fluid dynamics.
There have also been disclosed techniques for automatically generating a corresponding calculation program from a partial differential equation, an object region for the simulation and so on in C. Konno et al "Automatic Code Generation Method of DEQSOL", Journal of Information Processing, Vol 11, no. 1, 1987 and C. Konno et al "A High Level Programming Language for Numerical Simulation:DEQSOL" Denshi Tokyo (IEEE Tokyo Section), 25 (1986), in addition to the foregoing U.S. Pat. No. 4,841,479 (U.S. Ser. No. 361,283 and the document by Sagawa et al.
Those conventional program generating methods imply the following problems with respect to the fluid analysis:
(i) As a first problem, with the conventional program generating method, discrete points can be placed only on mesh points, which cannot deal with staggered meshes which are generally used for the fluid analysis problem. More specifically, for effecting a simulation of a fluid, variables are located at various points, however, since the conventional methods can only treat locations at mesh points, it is not possible to apply a method of placing discrete points at intermediate points between mesh points, whereby the solution is not converged.
(ii) As a second problem, the conventional program generating methods only prepare an ordinary central difference where discrete values are taken in the same proportion in the upstream side as well as the downstream side. Stated other way, the upwind difference generally used for fluid analysis problems has not been supported in the conventional generating method. A fine mesh, if can be formed, may attend to the upwind difference, however, the conventional methods cannot but form a rather rough mesh. Further, a generally performed discretization of the upwind difference is not regarded as a method for evaluating a value on the boundary of a control volume on the upwind side, and therefore, a method which replaces a differential operator as it is with a discrete equation of the upwind difference has been taken in general. In other words, the conventional methods do not take in consideration the foregoing upwind difference method, that is, a method for choosing a proportion of discrete values to be larger on the upstream side, and a central difference method has been utilized. As a result, it is not ensured that the input and output in a control volume in the vicinity of the boundary are identical. According to the theory of the control volume, an input quantity and an output quantity of a mesh point must be equal.
FIG. 23 is an explanatory diagram of an upwind difference method. If a third order upwind scheme according to Kawamura and Kuwahara is used in FIG. 23, by way of example, a first order upwind scheme, with a degraded accuracy, is used in a control volume 641 close to the boundary, while the third order upwind scheme according to Kawamura and Kuwahara is used only in an inner control volume 642. The discretization thereof is given by the following formulae: ##EQU1##
The formula (1) uses the first order upwind scheme while the formula (2) the third order upwind scheme, whereby physical input and output quantities on a control volume 643 are not equal to each other.
(iii) As a third problem, with the conventional program generating methods, the differential operators can be described only in previously prepared forms of .differential./.differential.x, .differential./.differential.y, .differential./.differential.z and a combination thereof. Also, since a central difference is employed as a method for discretizing these differential operators, the user cannot partially improve the discretization accuracy. More specifically, the conventional methods use discrete equations derived by the Taylor expansion, wherein the central difference is used in inner points while a one-side difference is used on the boundary. The central difference can provide a higher approximation accuracy, however, application of the central difference on the boundary results in that discrete reference points overflow outside of the region.
(iv) As a fourth problem, when a partial differential equation is in the form of an assignment statement such as "Unknown Quantity=Partial Differential Equation" solutions derived by the conventional methods are not always correct. A calculation program constructs a DO loop for mesh points effecting the same calculation when values at mesh points are to be calculated. If variables in the left side are also used in the right side, former values and new values may be mixed in points which are referred to when an assignment statement is discretized. Also, different orders of constructing the DO loops cause each reference point to use a different value, i.e., a former value or a new value, so that the same result cannot always be obtained.
(v) As a fifth problem, when a partial differential equation is solved by the use of a matrix solution, the conventional methods generate calculation programs which use matrix solutions for five-point difference in a two-dimensional space and for seven-point difference in a three-dimensional space, irrespective of the form of matrix, whereby it cannot be ensured that a solution appropriate to the form of a given matrix is always used.
(vi) As a sixth problem, the conventional program generating methods do not consider at all for the case where a physical boundary exists inside of an object region for the simulation and therefore cannot deal with such a problem.
(vii) As a seventh problem, the conventional program generating methods only-provides a method for generating a DO loop for a partial differential equation for performing a standard discretization and accordingly cannot deal with a partial differential equation which includes an upwind difference or a user-defined non-standard differential operator.
(viii) As an eighth problem, the conventional program generating methods effect a division by the dimension and volume of the control volume for each term in the right side of a partial differential equation after discretization, which results in increase in the number of operations and the length of equation, and mixture of errors due to the division.